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RUNAWAY STARS AND THE ESCAPE OF IONIZING RADIATION FROM HIGH-REDSHIFT GALAXIES 
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ABSTRACT 

Approximately 30% of all massive stars in the Galaxy are runaways with velocities exceeding 30 kms -1 . 
Their high speeds allow them to travel ~ 0. 1 — 1 kpc away from their birth place before they explode at the end 
of their several Myr lifetimes. At high redshift, when galaxies were much smaller than in the local universe, 
runaways could venture far from the dense inner regions of their host galaxies. From these large radii, and 
therefore low column densities, much of their ionizing radiation is able to escape into the intergalactic medium. 
Runaways may therefore significantly enhance the overall escape fraction of ionizing radiation, / esc , from small 
galaxies at high redshift. We present simple models of the high-redshift runaway population and its impact on 
/ esc as a function of halo mass, size, and redshift. We find that the inclusion of runaways enhances / esc by 
factors of ss 1.1 — 8, depending on halo mass, galaxy geometry, and the mechanism of runaway production, 
implying that runaways may contribute 50-90% of the total ionizing radiation escaping from high-redshift 
galaxies. Runaways may therefore play an important role in reionizing the universe. 
Subject headings: reionization — galaxies: high-redshift 



suggests that massive stars provide the bulk of the ionizing 
luminosity responsible for reionization (Mada u et ; 

m. S 



1. INTRODUCTION 

Observations of the cosmic micr owave background 
(iDunkley et al.l2009tlKomatsu et al|201 lh and the absorption 
spectra of high-redshift quasars dFan et al.l 12006) imply that 
the universe was reionized sometime between 6 < z < 12. The 
rapid drop in the quasar luminosity function at high redshift 

: io nizing 
dl 119991 

lLoeb & Barkanal l200Tt iBolton et al.ll2005l) . Stars at high red- 
shift are born within dense proto-galactic fragments, and thus 
much of their ionizing luminosity is thought to be attenuated 
by neutral hydrogen. The fraction of their ionizing luminos- 
ity that escapes into the intergalactic medium (IGM) and is 
available to reionize the universe, /esc? is a key quantity in 
determining how reionization occurred. 

Models for reionization that reproduce the observed UV 
luminosity function of high-redshift galaxies require both 
/esc ^ 0.1 a nd significant contributions from very faint 
galaxies (e.g.JPawlik et al . 2009; Srbino vskv & Wvithel2010| ; 
Haardt & Madau! 120121; iKuhlen & Faucher-Giguerd \20l% 
Bouwens et al. 201 laj). These requirements place very de- 
manding constraints on models of high-redshift galaxies. 

In the local universe, observations of both normal 
star forming and starburst galaxies favor escape frac- 
tions in the range of a few percen t (jLeitherer et al.l [T995j; 
Bland-Hawthorn & Malonevl Il999t iHeckman et alT I200U 
Grimes et al.l 120071) . The escape fraction seems to be sim- 
ilarly low at z ~ 1 (ISiana et al.l 120071: ICowie et all 120091: 
ISiana et aDl2010h . At higher redshifts (z ~ 3), observations 
suggest that / esc may increase to ~ 10%, at least for the 
relatively massive, highly star-forming galaxies for which 
measurements of f esf have been made dSteidel et al.l | 200~T| : 
Giallongo et alJl200a IShaplev et alJl2006h llnoue et aljbooa 
Siana et alJl2007t llwata et alJl2009t ISiana et alJl2010h . How- 



ble, owing to the high opacity of the IGM. Our understand- 
ing of how reionization proceeded therefore depends crucially 
on models and simulations that follow the escape of ionizing 
photons. 

Simple theoretical estimates of / esc can be obtained 
from analytic mod e ls of galaxies and HII regions 
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ever, constraints on / esc at z > 2 from gam ma-ray burst af - 
terglows favor low values of / esc ~ 0.02 dChen et alj|2007l) . 
While the data suggest that escape fractions may increase out 
to z ~ 3, it is worth noting that direct measurements of ion- 
izing photons during the epoch of reionization are impossi- 



(Dove&Shull 1994; Haiman & Loeb \199% iMadau et al 
1999t iDove et alj|2000t iRicotti & Shu !1[|2000t IWo od & Loeb 
2000; Fernandez & Shulll 1201 ll) . These models produce 
estimates of / esc that vary widely from 10~ 3 to ~ 1 depending 
on dark matter halo mass, dumpiness of the interstellar 
medium (ISM), gas density profile, and redshift. Sophis- 
ticated hydrodynamical sim ulations also predict a wide 
range of esc a pe fractions (Razoumov & Sommer-Larsen 
[20061 [20071; iGnedin etaT] 120081; IWise & Cenl 1200' 
Razoumov & Som mer-Larserj2010tlYaiima et al.ll201 ll) . 

While some of the differences between model predictions 
can be attributed to different choices for physical parameters 
such as halo mass and redshift, numerical methods and sub- 
grid models also contribute to the quoted discrepancies. For 
example, the geometry of the ISM and the relative distribution 
of stars and gas has a strong influence of / esc , and both are 
controlled by uncerta in prescription s for s tar formation and 
supernovae feedback. IClarke & Oeyl (120021) argued that there 
is a critical star formation rate above which the porosity of 
the ISM is approximately unity. A highly porous ISM is ex - 
pected to result in a high escape fraction dFuiita et al.ll2003l) . 
Indeed, simulations are able to achieve high escape fractions 
only when stella r feedback is ve ry strong and hence the poros- 
ity is very high. Ricotti (2002) has considered globular clus- 
ters as efficient sources of ionizing photons due to their high 
star formation efficiencies and high porosities. In this paper 
we consider an alternative scenario that is capable of produc- 
ing very high escape fractions without appealing to a highly 
porous ISM: the contribution from runaway stars. 

A significant fraction of massive stars in the Galaxy are 
moving at high velocity (> 30 kms" ) and are therefor e 
known as runaways (e.g.. lBlaauwfT96ltlGies & Boltonll 986). 
Two formation channels have been proposed for the formation 
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of runaways: 1 ) via dynamical ejections from young dense 
stella r systems ([Gies & B oltonll 1 9 8 6t iFui ii & Portegies Zwartl 
1201 U iPerets & Subri 120121) and 2) due to the explosion of 
a companion star (IZwickvl 119571: lBlaauwl|196lHGott etal.1 
Il970t iPortegies Zwartl 12000: lEldridge et all 1201 ll) . A com- 
bination of these two scenarios ha s also been considered 
(Pflamm-Altenburg & Kroupa 2010). Tracing the orbits of 
runaways back in time has provided eviden ce that both chan- 
nels c ontribute to the runaway population (Hoog erwerf et alj 
120011) . The available data on runaways allow for a wide range 
in the fraction of massive stars that are runaways, / mn , from 
0. 1-0.5 depending on definition an d correction for complete- 
ness and observational biases (e.g.. [Gies & Boltonll986tlGiesl 



fT987t ISt^JT99lt iTetzlaf f et aUkoill) . 

The goal of this paper is to assess the effect of runaways on 
the escape of ionizing radiation from high-redshift galaxies. 
Qualitatively we expect them to be important when galaxy 
sizes are smaller than ~ 100 pc, as then even relatively slow 
runaways with short main sequence lifetimes can travel from a 
galaxy's center to its outskirts before exploding. In this paper 
we will explore the effect of runaways by building simple ana- 
lytic models of high-redshift galaxies and modeling runaways 
formed through both the dynamical and supernova mecha- 
nisms. The potential importance of runawa y stars on the 
escape of ionizing photons was first noted by iDove & Shulll 
( 1994) in the context of the Galaxy. 

Throughout this paper we adopt cosmologi cal parameters 
consi stent with the 7th year WMAP estimates dKomatsu et alJ 
1201 ll) . namely = 0.73 and Sl m = 0.27, and a Hubble con- 
stant of H = 70 km s" 1 Mpc" 1 . 

2. MODEL 

Our goal is to construct simple models for the escape 
fraction of ionizing radiation from high-redshift galaxies 
both with and without the contribution of runaways. Our 
basic setup is motiv ated by and follows the approach of 
iRicotti & Shulll (120001) . and references therein. 

2.1. Galaxy Models 

We begin by considering the relation between cold dark 
matter halo mass and vi rial radius derived f rom cosmologi- 
cal A^-body simulations dNavarro et al.lll997l) : 
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(1) 

where fl m and f2 m (z) are the matter densities of the universe at 
the present epoch and redshift z, respectively, in units of the 
critical density (p c rit = 3Hq/8/wG), and A f is the overdensity 
thresh old of virialized dark ma tter halos, which depends on 
fi m (z) dBryan & Normanlll998l) . For our fiducial model we 
adopt z = 10. 

We will assume that 100% of the baryons are in the cold gas 
phas e (i.e., the stellar fraction is sa at thes e early epochs; see 
e.g.. IRicotti et alj|2002tlwise & Cenll2009l) . implying: 



M gas = /fcMhalo, 



(2) 



where fa = -17 is 
(Komatsu efal] 120111) . 



the universal baryon fraction 
The relation between the gas 



scale radius and the halo virial radius is taken to be dMo et al.1 
1998): 



r„ = = r v 



V2f b 



(3) 



where A is the dimensionless spin para meter, which for halo s 
in simulations has a mean value of 0.04 (BullocketaiJ[2001). 
The galaxy angular momentum is a fraction jj of that of the 
halo. We assume jd/fb = 1 and A = 0.05, as these parameters 
provide a good fit to the obs erved size distribution of local 
galactic disks dMo et al.ll 19981) . 

High-resolution cosmological hydrodynamic simulations 
produce a wide range of galaxy morphologies at high 
redshift, ranging f rom spherical to well-ordered disks 
dWise & Cenl l2009t IPawlik et all l20llt iRomano-Dfazetal] 
1201 U iWise et all 120121) . Because the true nature of these 
galaxies remains unknown, and at present unobservable, we 
consider two galaxy geometries, spherical and disk-like. In 
our spherical model the gas has an exponential number den- 
sity profile, n{r) = no e~ T ' rg . In our disk model, we consider 

a gas profile with n(r,z) = noe~ R l r *e~ z lh, where z g is the disk 
scale height. We consider models where the disk height-to- 
disk scale length ratio varies with mass as: 
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The mass dependence of z, I r g follows from th e assumption of 
an isothermal disk (e.g., Wood & Loeb 2000). To ensure that 
z g /r g never exceeds unity we sum its inverse in quadrature 
with one: (z g /r g T 2 -> (z g /r g )~ 2 + l. 

In our fiducial model massive stars and stellar clusters are 
born tracing the gas density profile. We will also consider 
a model where the stellar density is proportional to the gas 
density squared (p* oc p 2 ) as in IWood & Loebl {2000). The 
slope of the observed Kennicutt-Schmidt relation between star 
formation rate surface density a nd gas surface d ensity falls 
in between these two extremes (Kennicutt 1998). The gas 
and stellar density distributions are truncated at 5r g = 0.2r v ; r 
(and 5zg for disks). In summary, for our fiducial set of model 
parameters, the halo mass uniquely specifies the distribution 
of gas and stars within the halo. 

2.2. Photon Escape Fraction 

Our next task is to compute / BS1 .. As d iscussed in 
IDove & Shulll (fT99l and IRicotti feShul d20(p, the HII re- 
gions created by massive stars will not be spherical if the local 
gradient in the gas density profile is large. While the HII re- 
gions are allowed to be non-spherical in our model, the equa- 
tions presented in this section assume that the gas distribution 
is spherical; the generalization to disk geometries is straight- 
forward. Consider a star a distance r from the center of the 
galaxy. A point residing a distance r' from the star at an angle 
8 from the line connecting the galaxy center and the star will 
be at a distance R = \/r 2 + r' 2 + 2rr'cos9 from t he center of 
the galaxy (see Figure 1 in IRicotti & Shu ll 2000, for a sketch 
of the geometrical configuration^ . At each angle 8 one can 
compute the quantity: 



f(0) = l 



4-iraf] 



n H {R) 2 r l2 <\r', 



(5) 



where Q = Qh is the hydrogen ionizing luminosity, is 
the hydrogen case B recombination coefficient, and tin is 
the number density of hydrogen. We adopt = 2.6 x 

1 There appears to be an erroneous substitution of sin for cos in their equa- 
tion for R. 
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10~ 13 cm~ 3 s _1 , appropriate for gas at T = 10 4 K. Values of 
/ > indicate a non-zero escape fraction, while / < indi- 
cates that 100% of the ionizing radiation is absorbed. The 
critical angle for escape, 8 C is defined through the equation 
f(0 c ) = 0.0. Notice that 6 C is a function of Q and r. The angle- 
averaged fraction of ionizing photons escaping the galaxy is 
then computed by integrating f{8) over a solid angle: 



(/esc) = ^ / ' f(ff)2iramOM. 
4?r J 



(6) 



The total escape fraction of the galaxy is then computed by 
integrating over r and Q: 



/esc — 



(f e sc)P(r)P(Q)drdQ, 



(7) 



where P(r) and P(Q) are the probabilities of a star being at 
location r and with ionizing luminosity Q. 

We have neglected the effect of overlapping HII regions 
in our analytic framework. This is a common simplifica- 
tion made in the modeling of the ionizin g escape fraction 
dDove & Shull|[T994tlRicotti & Shulll2000b . The effect of in- 
cluding HII overlap would be to increase / esc , although the 
low e xpected star formation rates in z = 10 halos dWise & Cenl 
2009) suggests that HII overlap may not be significant. 

In the absence of runaways, P{r) and P(Q) are straightfor- 
ward to spec ify ( our prescription for including runaways is 
described in cj2.31 l. The former is determined by the gas den- 
sity profile, while the latter is taken to be the observed ion- 
izing luminosity function (LF) in local galaxies. The LFs 
of OB associations in nearby galaxies have power-law in- 
dices that range between —1.5 and -2.5 over the range 48 < 
log(Q AT 1 ) < 51 Jkennicutt et ail [19891 iMcKee & Wi lliams 



1997; Whitmore et al. 1999; Ivan Zed 120001: lAzimlu et all 
20111) . lOev & Clarke! ( 119981) argued that the data are con- 
sistent with a universal power-law slope of -2.0. We there- 
fore adopt an ionizing LF with a P(Q) oc Q~ 2 over the range 
48 < log{Q/s~ l ) < 5 1 . For reference, log(0 = 48 corresponds 
to a single 20 M Q star. We will explore variation in the lower 
luminosity cutoff, Q m i„, in S|3] Notice that our treatment of 
P(Q) naturally accounts for the fact that massive stars are clus- 
tered. 

Current evidence favors very little if any dust in low mass 
galaxies at high redshift, as in ferred from their very blue UV 
colors (Bo uwens et al.1201 lbh . The lack of dust is also consis- 
tent with their expected low metallicities. Our fiducial model 
is therefore dust-free, but we nevertheless explore the effect 
of dust on /esc- To incorporate dust into our models, we a dopt 
a wav elength depend ent dust cros s sec tion, a(\) from iPeil 
(119921) . as updated in lGnedin etail (12008b . Essentially noth- 
ing is known about the dust composition nor grain size dis- 
tribution in high-redshift galaxies, so we consider both Large 
and Small Magellanic Cloud-type dust (LMC and SMC, re- 
spectively) to give a sense of the variation expected from 
different grain populations. To make the computations sim- 
ple, we have taken averages of a(X) over ionizi ng photons, 
weighted by model stellar flux distributions from Smi th et all 
(2002). The average depends weakly on the effective tem- 
perature, varying by ~ 10% over the relevant range. The re- 
sulting flux-weighted mean cross sections for LMC and SMC 
dust are a = 1 x 10~ 21 cm -2 and a = 5 x 10~ 22 cm" 2 , respec- 
tively. For each r, 9 pair we can compute = aNn, where 
Nh is the column density of hydrogen (we adopt the observed 



dust-to-gas ratios of the LMC and SMC, which likely overes- 
timates the influence of dust at high redshift). The hydrogen 
ionizing luminosity Q that enters into Equation[5]is then sim- 
ply attenuated by e~ Td . 

With the gas density profile, ionizing LF, and dust model 
specified, / esc for non-runaway stars is then uniquely deter- 
mined as a function of halo mass. We now consider the addi- 
tion of runaways. 

2.3. Runaways 

For the purposes of computing / esc , runaways differ from 
non-runaways in two respects: they travel from their birth en- 
vironment, resulting in a change to P(r), and the runaway pop- 
ulation may not reflect the overall massive star population in 
its mass distribution and therefore its ionizing LF. As noted 
in the Introduction, there are two mechanisms capable of pro- 
ducing runaways in the Galaxy: 1) dynamical encounters in 
young dense stellar systems, and 2) the explosion of a close 
companion. In order to explore the effects of runaways on 
/esc we develop simple models for both scenarios, motivated 
by the observed runaway population in the local universe. 

In the dynamical mechanism, runaways inherit a high ve- 
locity, V mn , not long after their birth. We model runaways 
produced via this mechanism by taking a fraction of all mas- 
sive stars, /run, and giving them a velocity V ran at birth with 
a random direction in the galaxy. The runaway is allowed to 
travel for a time equal to its main sequence lifetime. The ef- 
fect of the potential due to the galaxy and halo mass on the 
motion of the runaways is neglectecfl 

In the simplest version of this model we assume that the 
runaways are, at birth, a random sampling of the overall mas- 
sive star population. However, in reality one may expect the 
maximum mass of runaways produced in this mechanism to 
be lower than the upper mass of all stars. This may occur 
because dynamical en counters tend t o eject the lowest mass 
star in the encounter (Anosova 1986). We model this effect 
by considering a possible truncation in stellar mass, M max , of 
the runaway mass function. 

We explore models where V mn either takes a sin- 
gle value or a distribution. T he numerical model s of 
iFujii & Portegies Zwartl(l201 lb and lPerets &"S ubr (|2012|) pro- 
duce runaways with a power-law distribution of V ran with 
an index of « — 1 and ps —1.5, respectively, over the range 
20 kms -1 < y mn < 100 kms -1 . Observations of runaways in 
the Galaxy seem to favor a Maxwellian velocity distribution 
with a dispersion of 30 kms" 1 dStondll991h iTetzlaff et all 
1201 ll) . implying a mean runaway velocity of ss 45 kms -1 . We 
consider both of these distributions and their effects on / esc . 

For the dynamical model we assume that the runaway LF is 
equal to the non-runaway LF up to a maximum ionizing lumi- 
nosity corresponding to a mass of 100 M . Conceptually, this 
corresponds t o replacing clusters w ith single stars. See the 
discussion in lOev & Clarkd ([1998), where they find that the 
single-star LF does not significantly differ from the cluster LF 
for a range of assumptions. To calculate the influence of run- 
aways, we must translate a given luminosity into a stellar mass 
so that we may assign stellar lifetimes. We have constructed 

2 The circular velocity of dark matter halos at 0.1r v j r « 3r g is < 30 kms~' 
for M hato < 10 9 M Q and < 60 kms" 1 for M halo < 10 10 Mq INavarro et all 
1997). For the typical runaway velocities we consider (Vmn > 30 kms - '), the 
effect of the halo potential on the runaway trajectory can be safely neglected 
at low halo masses, but it may have a modest effect at the highest masses we 
consider. 
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an ionizing luminosity-stellar mas s relation, Q{m), by com- 
bining the stellar interior models of Sch aller et al.l (11992b with 
the stellar spectral library of Smi th et all (120021 ) at 0.05Z Q . 
We then use t he main-seque n ce life time-mass relation, (m), 
presented in Ekstro m et al.l (1201 lb for their rotating stellar 
models. At fixed mass, the main sequence lifetime varies by 
±20% dep ending on the modeling of the stellar interiors and 
metal licity dMarigo et a0200lUSchaerel2002llEkstrom et alj 
l20ll . 

The second mechanism we consider is the supernova ori- 
gin of runaways. In this scenario the runaway was at some 
point in a binary with a more massive companion. When the 
companion explodes, the surviving star inherits some frac- 
tion of the pre-explosion orbital velocity. As in the previ- 
ous scenario, we specify the fraction of runaways, f mn , pro- 
duced in this way. In addition, we must specify the distri- 
bution of mass ratios, q = M2/M1, which we assume to be 
flat over the rang e 0.1 < q < 1.0, con sistent with constraints 
from the Galaxy (ISana & Evans 20TH) . We also assum e that 
the primary mass, M\ is drawn from a lSalpeterl dl955) initial 
mass function. Calculation of the runaway velocity requires 
knowledge of the period distribution, initial-final mass rela- 
tion for the more massive companion, and, for binaries that 
remain bound, kn owledge of the post-explosion orbital prop- 
erties (Gott[ |1972h . Consideration of t hese factors is beyond 
the scope of the pres ent work (see e.g., Portegies ZwartCOOOt 
lEldridge et al.l I20TTI) . Instead, we take an approach that is 
similar to our treatment of V ran in the dynamical model: we 
simply consider a fixed value of V mn that is similar to the ob- 
served average runaway velocity of Galactic O stars. 

The supernova mechanism differs from the dynamical 
mechanism in at least two important respects. First, the life- 
time of a runaway is shorter in the supernova scenario because 
the runaway is formed only after the death of a more massive 
companion, while in the dynamical scenario the velocity kick 
is imparted at birth. As we will see below, this effect can be 
significant. The ratio in median runaway lifetimes between 
the two mechanisms is ~ 3, so runaways produced via close 
dynamical encounters can have a substantially larger impact 
on /esc- In reality, runaways formed in the dynamical channel 
will be born some finite time after their birth, tempering this 
difference. The second distinction is that the runaway ioniz- 
ing LFs are different in the two models. While in the dynam- 
ical model runaways have the same LF as the non-runaways, 
in the supernova model the runaway LF is determined by a 
combination of the IMF, main sequence lifetime-stellar mass 
relation, and binary properties. In practice, the supernova run- 
away LF is not very different from the dynamical runaway LF 
— it is approximately a power-law with an index of ~ —1.7. 

Both of these models are implemented via Monte Carlo 
simulations. For our fiducial model we adopt the dynamical 
mechanism for runaways with f mn = 0.3, V mn = 60 km s _I , and 
M max = 100 Mq. In addition to this fiducial model, a variety 
of model permutations will be considered in $3] 

2.3.1 . Runaways at High Redshift 

We have adopted runaway models based on observations of 
massive stars in the Galaxy, so it is reasonable to ask whether 
the population may be different at high redshift. For both 
the dynamical and the supernova model, the binary popula- 
tion affects the runaway fraction, with more binaries leading 
to more runaways. High-mass binary formation may be an in- 
evitable cons equence of disk fr agmentation driven by high ac- 
cretion rates (Kratte ret al.ll2.01Qb . Indeed, recent simulations 



1.000 




FIG. 1. — Distribution of ionizing escape fractions for individual massive 
stars in the fiducial model for Mh a i = 10 s Mq. The non-runaways (black) 
and runaways (red) are shown separately. 

of population III star formation suggest that binaries may be 
ubiquitous even in zero metallicity systems (Turk et al1 l2009t 
Greif et al. 201 1). It is therefore plausible to assume that the 
binary fraction at high redshift is similar to that found in the 
local universe. 

In the dynamical ejection model the number of run- 
aways depends on star clu ster properties. Following 
Fuiii & Portegies Zwart (201 1), who suggest that clusters pro- 
duce a fixed number of runaways independent of cluster mass, 
one expects the runaway fraction to scale inversely with clus- 
ter mass. Therefore, if more stars are born in smaller —N clus- 
ters, a higher fraction will be runaways. Simple estimates 
of the characte ristic fragmentation mass, M c , in a Q = 1 disk 
( Toomre1 [l964l) suggest that it should be smaller at high red- 
shift. In a gaseous disk M c scales as M c ~ Sit 2 /k 2 where X 
is the mass surface density, a is the velocity dispersion in the 
gas, and n is the epicyclic frequency. Adopting Keplerian ro- 
tation and using Equations Q] and [2] yields M c ~ er 3 (l +z)~ 3 / 2 . 
Unless a is considerably higher in high-redshift disks com- 
pared to the Galaxy, we expect M c to be lower at high red- 
shift. Everything else being equal, this would suggest that the 
runaway fraction may increase with redshift. However, the 
efficiency of producing runaways must decline for the lowest 
cluster masses, where the number of high mass stars that can 
form is limited by the cluster mass. 

Given all these considerations, we believe that adopting a 
minimal model, where the runaway fraction does not evolve 
to high redshift, is acceptable. 

3. RESULTS 

Using the models described above we calculate both the dis- 
tribution of / esc values for massive stars in a given galaxy as 
well as the galaxy-averaged / esc for a range of halo masses. 
We begin with the results from our fiducial model, and then 
explore the dependencies of / esc on a variety of parameter 
choices. 
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FIG. 2. — Top Panel: Escape fraction of ionizing radiation, / eS c, as a func- 
tion of halo mass, Mh a | , for our fiducial set of model parameters. The escape 
fraction is shown for all stars, runaways only, and non-runaways only. Bot- 
tom Panel: The enhancement in / e sc , defined as the ratio between the escape 
fractions computed with all stars and with non-runaways only. Vertical lines 
represent the approxi mate halo masses of galaxies with an absolute UV mag- 
nitude of -12 and -16 (Kuhlen & Faucher-Giguere 2012). 

3.1. Fiducial Model 

We show the distribution of / esc values for stars in our fidu- 
cial model in Figure Q] for both runaways and non-runaways 
at a halo mass of 10 8 Mq. The vast majority of non-runaways 
have / esc = 0.0, with only a few per cent in the "transluc ent" 
regime where < / esc < 1 (see also Gnedin et al. 20QU). In 
contrast, a significant fraction of runaways have very high 
/e S c, with ss 50% having / esc > 0.6. No stars have / esc = 1.0. 
This is because stars are never infinitely far from the galaxy, 
so the galaxy always occupies a non-zero solid angle capable 
of absorbing some ionizing radiation. 

Figure [2] contains the main result of this paper. The top 
panel shows / esc as a function of Mhaio both for the galaxy as 
a whole and separately for the runaways and non-runaways. 
The escape fraction for non-runaways is < 3% for Mhaio > 
10 8 Mq, in agre ement with previou s analytic models of galax- 
ies at z = 10 (Wood & Loeb 2000). For non-runaways, / esc 
increases toward lower masses because the ratio of surface 
area to volume increases toward smaller galaxies (/ esc oc r~ l 



for non-runaways). The escape fraction for runaways, how- 
ever, is always larger than for non-runaways, and approaches 

1/3 

unity for the smallest galaxies (recall that r g oc M^ Q , making 
the distance traveled by runaways in low mass halos a larger 
fraction of the galaxy size). 

The bottom panel in Figure [2] shows the ratio between the 
overall escape fraction and the non-runaway value, which 
we call the enhancement factor in / esc . The turnover at 
Mhaio > 10 8 Mq is due to the fact that runaways travel on 
average a fixed physical distance that becomes a smaller frac- 
tion of the galaxy scale radius at larger masses. At sufficiently 
large scale radii the runaway escape fraction should converge 
to the non-runaway value, as the runaways will travel an in- 
significant fraction of a scale length. The runaway / esc must 
therefore decrease more rapidly than the non-runaway value. 
The enhancement factor has a maximum because / esc for run- 
aways begins to saturate near Mh a i = 10 8 Mq. At lower 
masses / esc for non-runaways continues to increase, so the ra- 
tio of /esc between all stars and non-runaways decreases at 
Mhaio <10 8 M Q . 

In this figure we also label approximate UV magni- 
tudes of galaxies in halos of masses 10 8 Mq and 10 9 M 
at z = 10 based on estimates of the z = 10 UV LF 
dKuhle n & Faucher-Giguej^ 120121) . These estimates are 
highly uncertain; they are included here merely to indicate 
the types of galaxies one might expect to occupy these halos 
at high redshift. 

The key result is that / esc increases by factors of 2 - 8 for 
10 7 Mq < Mhaio < 10 10 Mq when runaways are included in 
the model. Figure |2] demonstrates that, in our model, run- 
aways contribute 50-90% of all the ionizing photons that es- 
cape from the galaxy. Runaways may therefore play an im- 
portant role in reionizing the universe. 

3.2. Model Dependencies 

The sensitivity of these results to the assumptions of the 
model is explored in Figure [3] In this figure we vary the 
adopted redshift from 8 < z < 12, the dust attenuation, the 
minimum and maximum stellar mass, the adopted gas density 
profile, n(r), V ran , and the ionizing LF. We also explore the 
predictions of the supernova mechanism for the production of 
runaways. 

The trends are straightforward to interpret. At higher red- 
shift, galaxies are smaller at fixed mass. As the size of the 
galaxy decreases the effect of runaways becomes stronger 
both because it is easier for runaways to move to regions of 
low column density and because, at fixed galaxy mass, smaller 
sizes result in de nser galaxies, loweri ng the escape fraction for 
non-runaways (Woo d & Loebll20"o"ol) . For example, between 
z = 8 and z = 12 / esc increases by 70% at Mhaio = 10 9 Mq. 

Increasing V mn also results in runaways having a greater ef- 
fect on /-sc. In this figure we have also considered distribu- 
tions of V ran , rather than single values. We have included both 
a power-law distribution with an index of -1 and a mean of 
40 kms -1 , and a Maxwellian with a dispersion of 30 kms~ . 
In the latter model we have simultaneously set / mn = 0.46 in 
order to mimic the obser vationa l const raints on the runaway 
population determined by IStond (119911) . Adopting a distribu- 
tion of y ran rather than a fixed value has a modest effect on the 
results. 

Dust has a very minor effect on the derived / esc values. As 
pointed out by Gnedi n~et al.l (12008). even without dust the vast 
majority of non-runaways already have / esc = 0.0 (see Figure 
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FIG. 3 . — The effect of various model parameters and assumptions on the derived enhancement in / esc due to runaway stars. The solid line in each panel is our 
fiducial model where runaways are produced by dynamical encounters. Top Left: Variation in the redshift of the model from 8 < z < 12. Top Middle: Variation 
in the treatment of dust. Top Right: Variation in the maximum stellar masses of the runaways and in the minimum ionizing luminosity, Q m j„, of the combined 
runaway and non-runaway population. In the fiducial model, M max = 100 Mq and log(Q m j n ) = 48.0. Bottom Left: Variation in the gas density profile. The 
fiducial model of a spherical exponential profile is contrasted with an exponential disk with a Gaussian vertical density distribution. The disk model has a scale 
height-to-scale length ratio of 0.2 at Afhalo = 10 8 Mq and varies as M~ 2 / 3 . Bottom Middle: Variation in the treatment of the runaway velocities, Vmn- A fixed Vruo 
is varied from 40 — 80 kms -1 . In addition, we consider a power-law distribution of V rlin (labelled "POW" in the figure), and a Maxwellian distribution of Vmn 
as advocated in Stone 1 1991, labelled "S91" in the figure). In this last case / mn =0.46, as advocated by Stone. Bottom Right: Variation in the runaway model 
between our fiducial dynamical model and the supernova model. See the text for details. 



Q}, and so the addition of dust attenuation provides a minor 
modulation to the source-averaged / esc . In addition, the ma- 
jority of runaways reside in regions of very low column den- 
sity, where dust attenuation also has a minor effect. 

The adopted gas density profile has a dramatic effect on the 
resulting enhancement factor. In Figure [3] we show results 
for both a spherical exponential distribution and an exponen- 
tial disk with a Gaussian vertical distribution. Recall that in 
our disk model, the scale height-to-scale length ratio varies as 
M" 2 / 3 , as expected for an isothermal disk, with a value of 0.2 
at Mhaio = 10 8 Mq, In the case of disks, the runaway escape 
fractions are very high at high masses (/ esc > 50%) because 
the disks are thin and runaways can easily escape in the direc- 
tion perpendicular to the disk. At fixed mass, the enhancement 
in /esc is a strong function of Z g /r g , with smaller values result- 
ing in larger enhancements. At small masses the scale height- 
to-scale length ratio approaches unity in our model, and so the 
disks are essentially spheres. Thus the enhancement factor is 
similar between the two geometries^. 

Reducing the maximum stellar mass of runaways from 

3 The enhancement factors are not identical even at the lowest masses be- 



100 Mq to 60 M & causes a modest increase in the runaway 
enhancement because the mean stellar lifetime of runaways 
increases as the mass decreases. The runaways are therefore 
capable of traveling further from the central regions of the 
galaxy before they explode. We have assumed here that the 
runaway fraction does not depend on stellar massat> 6QMq. 
This assumption is supported by the results of IStond (119911) . 
who finds no variation in f mn amongst the O-type stars. When 
the minimum Q for both the runaways and non-runaways in- 
creases from 48.0 to 48.5, the decrease in the mean lifetime 
of the runaways results in a lower enhancement in / esc . 

The bottom right panel of Figure [2] shows the effect of 
adopting the supernova model for runaways rather than the 
dynamical mechanism. Supernova-produced runaways result 
in a very modest enhancement in / esc . The primary difference 
between these mechanisms, as implemented herein, is the dif- 
ferent lifetime of the runaways. Dynamically-created run- 
aways travel for roughly three times longer than supernova- 
created runaways, and so the former can more easily reach 
lower column densities in the galaxy. Observations favor 

cause a sphere is not identical to a disk with z g /r g fs 1 . 
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FIG. 4. — Variation in / C sc between all models considered in the text. The 
shaded regions encompass the maximum range of the models both with and 
without the inclusion of runaways (labeled "all" and "non-runaways" in the 
figure). Thick black lines represent the fiducial model. We have separately 
plotted the models that deviate substantially from the general trend, including 
models with a disk geometry, with supernova runaways, and with the stellar 
density profile scaling as the gas density squared. 

roughly equal contributions from dynamically-produced and 
supernova-produced runaways, and so we can expect reality 
to lie in between these two extremes. Various stellar evolu- 
tion models that predict different lifetimes for massive stars 
should change / esc by a more modest factor. 

There are additional parameters, not explicitly shown in 
Figure [2] that can also have a large impact on the enhance- 
ment in / esc . The runaway fraction, f mn , clearly will have a 
direct (linear) effect on / esc . The galaxy size also has a strong 
effect on the results, as can be inferred from the upper left 
panel of Figure [2] In this panel we consider variation in the 
adopted redshift, but this is simply changing the galaxy size at 
fixed mass (see Equation[T|i. The variation in redshift amounts 
to only a ±20% change in the typical size of a galaxy at fixed 
mass, implying that the results do depend sensitively on size. 
This highlights the need for incorporation of runaways into re- 
alistic simulations of high-redshift galaxy formation to assess 
their effect on / esc in a more quantitative manner. 

In Figure [4] we show / esc as a function of Mh a i for the full 
range of models discussed in this section. The shaded re- 
gions encompass the total variation in f esc due to the different 
model assumptions. Results are shown both with and without 
the inclusion of runaways. In the case where runaways are 
included, we separately highlight the most deviant models. 
These models are the disk geometry model and the supernova 
mechanism for the creation of runaways. In the disk model, 
/ esc is much higher than the fiducial model at high masses 
because runaways can more easily escape the galaxy when 
traveling perpendicular to the disk. In the supernova runaway 
model /esc is considerably smaller because runaways travel for 
a relatively short time before they explode. In all cases where 
runaways are produced via dynamical encounters the escape 



fraction can be quite high even in moderately large halos. 

Finally, in Figure |4] we separately show a model where the 
stellar distribution is proportional to the gas density squared, 
in contrast to our fiducial model where the stellar density is 
linearly proportional to the gas density. This model takes 
into account that stars form preferentially at the highest den- 
sities. In this model the escape fraction for non-runaways is 
10-20 times lower than our fiducial model, as the stars are 
now much more embedded on average. On the contrary, the 
escape fraction of runaways is identical to the fiducial model 
for Mhaic. ^ 10 9 Mq. In larger halos the galaxy size becomes 
comparable to, and eventually larger than, the mean distance 
traveled by runaways, and so the birth environment of the 
runaways becomes increasingly important in determining / esc . 
Since the non-runaway / esc is so much smaller than the run- 
away / esc in this model, the behavior of the runaways entirely 
determines the behavior of the overall escape fraction. As is 
clear from Figure |4] this model produces very similar overall 
/esc values for M < 10 9 M compared to the fiducial model. 
At higher masses this model produces lower escape fractions 
because the runaways are closer to their birth environment (in 
units of scale lengths), and bom in denser regions on average 
when oc p 2 g . 

4. DISCUSSION & CONCLUSIONS 

In this paper we have shown that the ionizing radiation from 
runaway stars may contribute substantially to the reionization 
of the universe. These stars migrate toward the low-density 
outer regions of high-redshift galaxies where their radiation 
can easily escape into the IGM. The importance of their mi- 
gration is enhanced at high redshift because the galaxies are 
much smaller than at z = (by a factor of ~ (1 +z) _1 ). Assum- 
ing that runaways have a prevalence that is similar to what is 
observed in the Galaxy (~ 30% for massive stars) and that 
dynamically-created runaways constitute a significant frac- 
tion of all runaways, they can increase the total escape fraction 
of ionizing photons from high-redshift galaxies by factors of 
2-8 compared to the escape fraction of non-runaway stars. 

Our conclusions depend strongly on three model ingredi- 
ents: the size of the galaxy, the runaway fraction, / mn , and 
the production mechanism for runaways. The effect of run- 
aways on /sc depends linearly on / ran , and therefore if the 
fraction of runaways is substantially smaller than what we 
have assumed here, their importance in an extragalactic con- 
text will be limited. We have also assumed that galaxy sizes 
at high redshift scale with halo size in a manner similar to 
what is found at z = 0, and so high-redshift galaxies are as- 
sumed to be much smaller than local galaxies. Small galaxies 
strongly enhance the effect of runaways on / esc - Finally, we 
have shown that runaways produced via dynamical encoun- 
ters have a much larger effect on / esc than runaways produced 
via an explosion of a close companion because of the differ- 
ent runaway lifetimes implied by these models. Observations 
of runa ways in the Galaxy fa vor a mixture of these two mech- 
anisms (Tetzlaf f et alj|201 H) . and so we can expect that run- 
aways in high-redshift galaxies will play an important role in 
the escape of ionizing radiation. 

The relative importance of runaways depends on the es- 
cape fraction of non-runaway stars. In our model the 
non-runaways have / esc < 10%, and so runaways can be 
very influential. However, if the non-runaway escape frac- 
tion were much higher, then the runaways would neces- 
sarily have a diminished impact. As mentioned in the 
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Introduction, some recent hydrodynamic simulations have 
fou nd very high escape fractions in high-redshift galax- 
ies dWise & Cenl 120091: iRazoumov & Sommer-LarseiilBOlOt 
Yaii ma et al.ll2011l) . These models invoke very efficient feed- 
back that leads to a highly porous ISM which in turn allows 
for many unobscured sight-lines to massive stars (see also 
IWood & Loebll2000t iFernandez & Shullll201 ll) . If these nu- 
merical models are correct, then the influence of runaways on 
/esc will be significantly decreased. To better understand the 
influence of runaways on realistic galaxies, they must be in- 
cluded self-consistently in detailed hydrodynamic simulations 
that include the effects of radiative transfer. 

The influence of runaways on / esc can also be tested directly 
with observations of high-redshift galaxies. If the escape frac- 
tion of non-runaways is indeed low, then there will be signifi- 
cant luminosity from recombination emission lines associated 
with the gaseous disk. Beyond a few disk scale lengths the 
runaways will dominate the stellar emission, which implies 
a decrease in the ratio of recombination line flux to ultravi- 
olet flux. In other words, the scale length of a galaxy mea- 
sured in recombination lines should be smaller than the scale 
length measured in broadband ultraviolet flux. Such an obser- 
vation would provide evidence for the influence of runaways 
(although we note that other explanations for a varying re- 
combin ation line-to-ultr aviolet flux ratio have been proposed, 
see e.g jLee et al.ll2009l) . The ratio of scale lengths of the run- 
aways to non-runaways is a decreasing function of halo mass: 
at Mhaio = 10 8 Mq the ratio of scale lengths is 2.6 while at 
MhMo = 10 9 Mq the ratio is 1.5, for our fiducial model. Such 
measurements should be within reach of the next generation 



of thirty-meter telescopes. In addition, the occurrence of su- 
pernovae far from the inner regions of high-redshift galaxies 
should be common, and would provide another diagnostic of 
the prevalence of runaways during the reionization epoch. 

Runaways may have other i mportant effects on the ev olu- 
tion of galaxies and the IGM. ICeverino & Klypinl j2009h in- 
cluded runaways in a high-resolution simulation of a Milky 
Way-like disk galaxy and found that they lead to much more 
efficient heating of the ISM. At high redshift, runaways will 
also be able to effectively heat and enrich the outer regions 
of halos and the IGM. Binary formation among the first stars 
may also have an important effect on IGM heating and chemi- 
cal enrichment via the production of Population III runaways. 

In this paper we have offered a plausible mechanism for the 
efficient escape of ionizing radiation from galaxies during the 
reionization epoch. Runaways almost certainly exist at high 
redshift, and their presence will result in an increase in / esc . 
Detailed simulations that incorporate runaways are required in 
order to make more quantitative statements regarding their po- 
tentially fundamental role in the reionization of the universe. 
Our results also highlight the need for more detailed observa- 
tions of the runaway population in the local universe, includ- 
ing the dependence of / ran on stellar mass, the runaway veloc- 
ity distribution, and the relative contribution of dynamically- 
produced and supernova-produced runaways. Fortunately, a 
wealth of new data o n runaway stars should be available in the 
coming years (e.g. JOev & Lamb1l201 UlLamb et al.ll201 ll) . 

We thank Andrey Kravtsov, Avi Loeb, and Sally Oey for 
comments on an earlier draft. We also thank the referee for a 
careful review of this manuscript. 
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